%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
load('counter.mat')

rng(1)

N_sim_div = 5000;

national_shares_sim = zeros(1,S);
for s=1:S
    national_shares_sim(1,s) = (sum(data_pat(:,s,end)) / sum(sum(data_pat(:,:,end))));
end

data_specialization = zeros(N,1);
for i=1:N

    data_specialization(i) = sum(((data_pat(i,:,end)/sum(data_pat(i,:,end))) - national_shares_sim).^2);
    
end

national_shares_sim = zeros(1,S);
for s=1:S
    national_shares_sim(1,s) = (sum(L_y_dyn(:,s,end)) / sum(sum(L_y_dyn(:,:,end))));
end

model_specialization = zeros(N,1);
for i=1:N

    model_specialization(i) = sum(((L_y_dyn(i,:,end)/sum(L_y_dyn(i,:,end))) - national_shares_sim).^2);
    
end

local_reliance = zeros(N,S,S);

for i=1:N  
    for r=1:S % Destination sector
        denom = 0;
        for s=1:S
            denom = denom + sum(lambda_dyn(:,s,end-1).*diffusion_frictions_dyn((s-1)*N+1:s*N,(r-1)*N+i,end).*(alpha_dyn(s,end)^theta));
        end
        
        for s=1:S % Origin sector
            local_reliance(i,s,r) = sum(lambda_dyn(:,s,end-1).*diffusion_frictions_dyn((s-1)*N+1:s*N,(r-1)*N+i,end).*(alpha_dyn(s,end)^theta)) / denom;
        end
    end
end

national_reliance = zeros(S,S);

for s=1:S
    for r=1:S
        national_reliance(s,r) = sum(L_y_dyn(:,r,end).*squeeze(local_reliance(:,s,r))) / sum(L_y_dyn(:,r,end));
    end
end

model_spec_adj = zeros(N,1);

for i=1:N
    for s=1:S % Origin sector
        term = 0;
        for r=1:S % Destination sector
            term = term + ((L_y_dyn(i,r,end) / sum(L_y_dyn(i,:,end)))*local_reliance(i,s,r) - (sum(L_y_dyn(:,r,end)) / sum(sum(L_y_dyn(:,:,end))))*national_reliance(s,r));
        end
        model_spec_adj(i) = model_spec_adj(i) + term^2;
    end
end

T_sim = 2;

rand_shocks = randn(S,T_sim,N_sim_div);

%%% Main model:

pop_sim = zeros(N,T_sim+1,N_sim_div);
pop_growth_sim = zeros(N,N_sim_div);

intra_id = 0;

for n_sim = 1:N_sim_div
    
    n_sim
    
    pop_sim(:,1,n_sim) = pop_dyn(:,end);

    alpha_sim = zeros(S,T_sim+1);
    alpha_sim(:,1) = alpha_dyn(:,end);
    
    perm_alpha_shocks = std(alpha_hat_dyn(:))*squeeze(rand_shocks(:,:,n_sim));
    
    lambda_sim = squeeze(lambda_dyn(:,:,end));
    L_k_sim = L_k_dyn(:,end); L_y_sim = squeeze(L_y_dyn(:,:,end));
    
    migration_exposure_sim = migration_exposure_dyn(:,:,end);

    for t = 1+1:1+T_sim
        
        L_k_sim_aux = L_k_sim;
        
        alpha_sim(:,t) = alpha_sim(:,t-1).*(1+perm_alpha_shocks(:,t-1));
        
        [pop_sim(:,t,n_sim),lambda_sim,L_k_sim,L_y_sim,L_o_sim,pi_sim] = counter_execute(N,S,d_geo,d_know,d_geo_distance,d_migr_const_aux,migration_exposure_sim,intra_id,lambda_sim,alpha_sim(:,t),squeeze(epsilon_dyn(:,:,end)),...
            pop_sim(:,t-1,n_sim),L_k_sim,L_y_sim,ubar_dyn(:,end),mu,theta,zeta,omega,fertility_dyn(end));
        
        prob_orig_given_dest_y_sim = zeros(N,N);
        for i=1:N
            for j=1:N
                prob_orig_given_dest_y_sim(i,j) = sum(pi_sim(i,j,:))*L_k_sim_aux(i) / sum(L_y_sim(j,:));
            end
        end
        migration_exposure_sim = prob_orig_given_dest_y_sim.*(1-eye(N)) ./ sum(prob_orig_given_dest_y_sim.*(1-eye(N)),1);
        
    end
    
    pop_growth_sim(:,n_sim) = log(pop_sim(:,end,n_sim)) - log(pop_sim(:,1,n_sim));
    
end

sd_pop_growth_sim = std(pop_growth_sim,[],2); 

for_stata = [cz_list, data_specialization, model_specialization, model_spec_adj, sd_pop_growth_sim, pop_dyn(:,end)];
csvwrite('results_model\spec_sim_for_stata_onebgp_temp.csv',for_stata)

%%% Intra model:

pop_intra_sim = zeros(N,T_sim+1,N_sim_div);
pop_growth_intra_sim = zeros(N,N_sim_div);

intra_id = 1;

for n_sim = 1:N_sim_div
    
    n_sim
    
    pop_intra_sim(:,1,n_sim) = pop_dyn(:,end);

    alpha_intra_sim = zeros(S,T_sim);
    alpha_intra_sim(:,1) = alpha_intra_dyn(:,end);
    
    perm_alpha_shocks = std(alpha_hat_intra_dyn(:))*squeeze(rand_shocks(:,:,n_sim));
    
    lambda_sim = squeeze(lambda_dyn(:,:,end));
    L_k_sim = L_k_dyn(:,end); L_y_sim = squeeze(L_y_dyn(:,:,end));
    
    migration_exposure_sim = migration_exposure_dyn(:,:,end);

    for t = 1+1:1+T_sim
        
        L_k_sim_aux = L_k_sim;

        alpha_intra_sim(:,t) = alpha_intra_sim(:,t-1).*(1+perm_alpha_shocks(:,t-1));
        
        [pop_intra_sim(:,t,n_sim),lambda_sim,L_k_sim,L_y_sim,L_o_sim] = counter_execute(N,S,d_geo,d_know,d_geo_distance,d_migr_const_aux,migration_exposure_sim,intra_id,lambda_sim,alpha_intra_sim(:,t),squeeze(epsilon_intra_dyn(:,:,end)),...
            pop_intra_sim(:,t-1,n_sim),L_k_sim,L_y_sim,ubar_dyn(:,end),mu,theta,zeta,omega,fertility_dyn(end));
        
        prob_orig_given_dest_y_sim = zeros(N,N);
        for i=1:N
            for j=1:N
                prob_orig_given_dest_y_sim(i,j) = sum(pi_sim(i,j,:))*L_k_sim_aux(i) / sum(L_y_sim(j,:));
            end
        end
        migration_exposure_sim = prob_orig_given_dest_y_sim.*(1-eye(N)) ./ sum(prob_orig_given_dest_y_sim.*(1-eye(N)),1);
        
    end
    
    pop_growth_intra_sim(:,n_sim) = log(pop_intra_sim(:,end,n_sim)) - log(pop_intra_sim(:,1,n_sim));
    
end

sd_pop_growth_intra_sim = std(pop_growth_intra_sim,[],2); 

for_stata = [cz_list, data_specialization, model_specialization, model_spec_adj, sd_pop_growth_sim,sd_pop_growth_intra_sim, pop_dyn(:,end)];
csvwrite('results_model\spec_sim_for_stata_onebgp.csv',for_stata)